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We study the interaction energy of the confining strings in the static rectangular tetraquark 
system in SU(3) gluodynamics. Using a two-state approximation we calculate the energy eigenvalues 
and corresponding eigenvectors of string states of different geometry. The interactions are studied 
both for the co-aligned and for the counter-aligned parallel strings as the functions of positions of 
the tetraquark constituents. We formulate a simple soap-film model and show that with a good 
accuracy the string interaction energy corresponding to the ground state of the tetraquark system 
can be described by this soap-film model. 
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I. INTRODUCTION AND MOTIVATION 
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, The complicated phenomenology of the strong interactions at low energies - such as the spectrum of mesons and 
baryons as well as the properties of their resonances - is determined by the forces between hadronic constituents, 
, quarks and antiquarks. These forces are dependent on spin orientations, relative velocities and distances between the 
^ > quarks and (anti-)quarks in the bound states P, Q- The major role in the QCD spectrum is played by the quark 
confinement which implies that the quarks are not observed as asymptotic states. 
■ In SU(3) gluodynamics the confinement manifests itself as a strong IGeV/fm) attractive force between an 
QQ ' external quark and an anti-quark as they get separated by large enough distances. It is believed that this force is a 
, result of a formation of a flux tube, or a QCD string, between the external color charges. Since the QCD string has 

a constairt energy per unit length the confining force is asymptotically separation- independent. 
\ In the real QCD the confining force is screened by pairs of quarks and anti-quarks which emerge from the vacuum 
and break an exceedingly long QCD string. Still the confining force caused by the QCD strings plays a crucial role 
in the structure of hadronic bound states and their resonances when these systems have sufficiently large spatial size. 
It is very plausible that the hadron structure is influenced not only by the interaction between the quarks themselves 



I 

^ ' {i.e., between the ends of the QCD string) but also by the interaction between the segments of the QCD string(s) with 
(~| , each other. The role of the QCD string interactions is the main subject of this article. Since the origin of the string 
and, consequently, the essence of confinement, is encoded in the gluonic degrees of freedom, we restrict ourselves to 
the quenched QCD case. 

■ The three-quark (baryon) state represents a simplest non-trivial case of the inter-string interactions. In the baryon 
5h I there are two possible geometries of the QCD string connecting three quarks: the A-like geometry which connects 
the quarks pairwise and the y-like geometry which connects quarks by three segments of the string with a junction 
in a middle. Lattice study of the baryon static potential revealed that the potential agrees with the dominance of the 
two-body interaction at small interquark separations |^ ^ |^ . This result is natural from the point of view point of 
the asymptotic freedom which predicts that at short distances the interquark interactions are reduced to the two-body 
interquark forces. 

The y-shape of the QCD string is difficult to observe at short distances due to the dominance of the Coulomb 
potentials ^. However, as the interquark separations get larger, the F-ansatz behavior prevails the two-body interac- 
tions null. The y-type profile of the QCD string was also observed at large inter-quark separations in quenched 
and unquenched lattice QCD in Refs. 0,0 where Abelian projection was used as a useful tool. In Ref. i8|| the F-type 
profile was observed in unquenched lattice QCD with the help of gauge invariant operators. 

The F-geometry of the fiux tube profile in baryon-like configuration of quarks is consistent with the prediction ^ of 
the dual superconductivity scenario 10], which is assumed to work in the infrared region (for a review see Refs. 
The field correlator method also predicts the F-geometry of the strings in the baryon p^ . 



^ Sec, however, Ref. |3, where an evidence of the y-profile at short distances was found. 
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In order to investigate the interaction between QCD strings it is convenient to consider a geometry with two 
long disconnected segments of the string. This corresponds to a generahzation of the baryonic system to a system 
consisting of four or more (anti)quarks. In this case, the interactions between the QCD strings can be deduced from 
the interquark potentials. The multi-quark potentials were intensively studied on the lattice "ff, Tl, El El E 113 • In 
the simplest case of the SU(2) gauge theory, the chromoelectric string docs not have a directional dependence since 
quarks and anti-quarks in this theory are indistinguishable from each other. A numerical study of energies in the 
four-quark systems shown that the interaction, or binding, energy is nonzero in the general case jl3 ]. It was found that 
if the four quarks are partitioned into a colorless two-meson state with basic states degenerate (or, almost degenerate) 
in energy, then the interaction energy of the system is strongly enhanced. It was thus shown that the energies of the 
four-quark states cannot be simply represented as a sum in terms of two-quark potentials. These features indicate 
that the mutual interaction of the chromoelectric strings can not be ignored in the SU(2) four-quark systems. 

In the SU(3) gluodynamics investigation of the interaction energy of the tetraquark system with two co-aligned 
mesons has shown (14. .15] that the strings form a joint profile with two Y-tyj>e junctions for small enough separations 
between the mesons. 

The two-meson interactions were also studied analytically in the dual superconductor model in the SU(2), Ref. 
and SU(3), Ref. HE! cases, as well as in the gauge-invariant field-strength cumulant approach 2, 20,]. 

Apart from purely theoretical interest in investigation of the inter-string interactions (potentially applicable to 
physical systems) there is a direct experimental motivation to study the 4Q systems in particular. Presently, candidates 
for a 4(3 state include the X{3872) resonance j2]J as discussed in Ref. 22], and Ds{2317) state 23]. The knowledge 
of the flux tube structure in the four quark systems seems to be important also in hadron scattering reactions since 
these reactions may also involve a rearrangement of the confining strings. 

In this paper we study the interactions of the QCD strings in the system consisting of two quarks and two anti- 
quarks. The general considerations concerning the interactions of the two (segments of) confining strings are discussed 
in Section m In order to calculate the interaction energy we diagonalize the QCD Hamiltonian in the basis of the two 
lowest energy states. As a guiding example, the calculation of matrix elements of a Hamiltonian is done analytically in 
Section Iml in a soap- film approximation which resembles the strong coupling expansion in the lattice gluodynamics. 
This calculation - which is of a qualitative nature - is then compared with the results of the lattice simulations of 
SU(3) gluodynamics discussed in Section Hvl Our conclusions are summarized in Section The technical details of 
the numerical procedures as well as principal output of simulations are presented in Appendices. 



II. INTERACTIONS BETWEEN QCD STRINGS 



The QCD strings have finite width but if they are long enough they can be considered as line-like objects connecting 
quarks or/and anti-quarks in three spatial dimensions. In a meson system the strings carry the chromoelectric flux 
which originates at a quark and is absorbed at an anti-quark. The interaction between the strings can be studied 
numerically with the help of the correlator of two large Wilson loops separated by a large enough distance. The 
correlator must be appropriately normalized by the expectation values of the individual Wilson loops, in order to 
subtract the energies of the infinitely separated quark-antiquark pairs^. The interaction energy Eint{d,r) is then 
defined as 

Eint{d, r) =- hm - log — ttt : 1 

t^oo t {Wrxt(0)) {Wrxt{d)) 

where Wrxt is the rectangular Wilson loop of size r x t, d is the distance between two parallel meson systems. The 
Wilson loop Wrxt describes the evolution of the quark-antiquark pair separated by the distance r from time = 
to time ti = t. 

In Eq. ^ the time extension t of the Wilson loops should be taken infinite. In order to study string interaction 
contribution into the full interaction energy one should also take r to be large in order to avoid effects of the short 
range Coulomb interaction between static quarks. In the real lattice simulations the above requirements are difficult 
to comply with. 



^ Note that flipping the orientation of a pair, r 



—r corresponds to the change W — > W t of the corresponding Wilson loop. 
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A. Flip-flop picture 



The important difference between the "ideal" strings discussed above and the "real" strings in lattice simulations 
is the fact that a lattice string has a finite length r and is evolved over a finite Euclidean time t. The constraints on r 
and t are due to the finiteness of the lattice volume. As a result of this finiteness, there are some effects which make 
the interpretation and realization of the Eq. Q difficult. 

The finite t implies the infiuence of higher excitations (this effect is common to all lattice studies) . We will address 
the question of the ^-dependence of our numerical results in Appendix IB II where we show that the results are rather 
stable against this kind of systematic errors. 

If the distance d is not large enough then the (anti)quark of one quark-antiquark pair is interacting with the one of 
the other quark-antiquark pair ^ providing an unwanted contribution to the interaction energy Ei^t- As a smearing 
technique we are using the hypercubic blocking (HYP) to improve the signal to noise ratio for our observables, also 
decreases effects of the short range inter quark interaction. As it was shown in Ref . |24^ , the perturbative interaction 
- which dominates the interactions of quarks at small separations - is suppressed by the HYP blocking. An interested 
reader can find the details of the HYP blocking in Appendix lA 21 

The effect of finite length r is seen in a string rearrangement which is also known as the " flip-fiop" effect discussed 
in the context of SU(2), Ref. and SU(3), Ref. 0|, lattice gauge theories. The fiip-flop is a string rearrangement 
process which happens as one changes the positions of the ends of the strings (quarks). The rearrangement is caused 
by a requirement for the strings to be in the ground state, i.e. in a state with lowest possible (for given boundary 
conditions) value of energy. As we will see later this is basically the same as the requirement for the string configuration 
to have a shortest possible length. 
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FIG. 1: A schematic illustration of the flip-flop picture: (a) the longitudinal and (b) the transverse string configurations. 



The rearrangement is illustrated in Fig. ^ for simplest case of two parallel strings with the same length and 
opposite orientations. This string configuration is called below the "longitudinal" configuration. The two quarks and 
two antiquarks shown in Fig. ^a) can also be connected by strings in another way shown schematically in Fig. ^b) . 
Below we call the latter "transverse" configuration. Let us keep the size r of the pair fixed and decrease the separation 
d. When d becomes smaller than some critical value dc the energy of the transverse configuration becomes smaller than 
the one for the longitudinal configuration. It is energetically favorable for the strings to realign at this point. Having 
neglected an interaction between the segments of the string, we can estimate an energy of the string configuration as 
E « crL, where L is a total length of the string configuration, and cr is the universal string tension'*. Obviously, due 
to a simple geometry of the string configuration shown in Fig. ^ the distances dc and r are approximately equal in 
this particular example. However, in other systems the competitive configurations can be more complicated, like the 
one shown in Fig. |21 for the case of equally oriented strings. The critical separation in this case is estimated to be 
dc ~ r/V3. 

If one treats the string configuration classically, the string re-alignment indicates existence at the fiip-fiop point 
d = dc of (at least) two degenerate energy minima in the space of the string configurations. On the quantum level, the 
degeneracy leads to the level splitting effect, as it is well known from the quantum mechanics. The level splitting is 
caused by the tunneling transitions between the wells of potential energy. Therefore one expects that the "transverse" 
and the "longitudinal" string states are mixed with each other due to the tunneling transitions and, consequently, the 
degeneracy is taken away. Technically, the tunneling transitions lead to an appearance of the non-diagonal matrix 
elements of the Hamiltonian in the basis of the original states (i.e., the states corresponding to the far-separated 



^ Note that the interactions of the quark constituents from the same pair are aheady excluded by the proper normaUzation of Eq. ITV . 
* The universahty of the string tension means that the string tensions in the heavy mesons and baryons are the same. This universaUty 
is supported by the numerical evidence obtained in lattice studies of the static baryon 0. 
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FIG. 2: The string configurations for tfie quark pairs witli tfie same orientation: (a) tfie longitudinal configuration; (b) The 
transverse configuration with two string junctions. The minimal length corresponding to the case (b) is Lj, = r + \/3d. 

"unperturbed" individual wells). To get the correct energy states one needs to diagonalize the Hamiltonian taking 
into account the tunneling transitions. The true ground state of the string in a tetraquark is a weighted mix of all 
"unperturbed" states of this system. In the two- and three-quark systems a similar investigation was performed, 
respectively, in Refs. 0| and [2^ in order to investigate gluonic excitations. 

In order to treat the degeneracy of the ground state exactly one needs to consider not only the ground state itself 
but also one should know all the states of the system. However, the higher is the state the lower is its mixing with the 
ground state. Therefore, one usually considers the two-state approximation, which takes into account only the ground 
state and the next energy state of the string system. Technically, this approach corresponds to the diagonalization of 
the Hamiltonian in the two dimensional subspace of the states. 

B. Two-level Approximation 

We can construct (" prepare" ) an operator to study the particular configuration of the quarks and anti-quarks by 
connecting the corresponding sets of points in a three-dimensional time-slice by gauge transporters. These trans- 
porters, called sometimes as "Schwinger lines", are the open Wilson lines in the fundamental representation of the 
gauge group. Such a set of lines transforms under gauge transformations as a set of quarks and anti-quarks. The 
better is the overlap of this operator with the real multiquark ground state (determined by the QCD dynamics) the 
smaller are effects of higher excitations. In order to improve this overlap we use the APE smeared gauge transporters, 
see Appendix I A II for details. 

We consider the configurations of the Wilson lines shown in Figures ^ and El The corresponding states define a 
two dimensional subspace of the state space discussed in the previous Subsection. The explicit definition of the basis 
of the states depends on positions of quarks and anti-quarks. We consider two simple geometries: co-aligned (Fig|21) 
and oppositely aligned (Fig^ parallel strings of the same length, and we present the exact expressions concurrently 
where the relative orientation makes a difference. 

The initial states of our basis - corresponding to the classical string configurations ~ are shown in Table The 
oriented lines mean the gauge transporters, the closed contour is the Wilson loop, and iVc = 3 is the number of 
colors used to normalize the initial states^. We refer to these basis states as |^) and \B) regardless of the orientation. 
Although the states are normalized, = {B\B) = 1, they are not, in general, orthogonal to each other, 

g^{A\B)^Q. (2) 

The basis states overlap g is proportional to the spatial Wilson loop and also shown in Table HJ 



^ Note that the number of the string segments coming to the string junction is equal to the number of colors. 
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TABLE I: The graphical definition of the states from the initial basis and of their overlap. 

Let us consider the diagonalization of the Hamiltonian H in the space |A)©|i3). To this end one can first diagonalize 
the Euclidean time evolution operator 



S{t) = exp(-fft) 



(3) 



the matrix elements of which can be calculated in lattice simulations. Despite the operator identity © is not correct 
in the finite basis, Sab ^ &^pi—Habt), one expects that only the lowest energy states contribute to the evolution 
operator matrix elements for large enough evolution time t. We will see that the t-dependence of the energy levels, 
defined via evolution operator eigenstates, has a form of l/t corrections and, consequently, can be neglected for large 
enough t. 

The main observables calculated in our numerical simulations are the S matrix elements: 



{A\e-"'\A) {B\e-^'\A) 
^ ' \ {A\e-"'\B) {B\e-"'\B) 



(4) 



They are expressed as the Wilson loops correlators which for the case of the oppositely aligned strings can graphically 
be represented as follows: 



Sit) = (0| 



V 

For co-aligned strings the graphical representation is a bit more complicated: 



|0> 



(5) 



S{t) = (0| 



N' 





|0) 



(6) 



\ ^ 1/ ^ / 

The problem to find eigenfunctions and eigenvalues of the matrix S 

S{t)\E) = eM~Et)\E) , 
can be reformulated as a generalized eigenproblem in terms of the matrix elements 

Sab{t)hb = X{t)Tabhb , 



(7) 
(8) 
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where a,b £ {A,B} and \E) — HaIA) + hB\B). The matrix T in the right-hand side of Eq. is nothing but the 
overlap of the basis states, discussed above, 

Equation Q together with the normahzation condition 

{E\E) =h\ + hl + 2ghAhB = l (10) 

defines two eigenstates |i?o,i) of the evolution operator (and, with all precautions, of the Hamiltonian) . These 
eigenstates correspond to the energies 

Eo,--'-^. (11) 

In order to solve the eigenproblem one should know all the matrix elements Satit)- To understand the qualitative 
properties of the eigenvalues we perform an analytical calculation in the soap-film approximation in the next Section. 

Note that we call the energy Eq as the ground state energy while the energy Ei is refereed to as the excited energy. 
This standard terminology - the ground state corresponds to the lowest energy while the excited state to the higher 
one - should not be taken literally. Indeed, the states correspond to different string configurations in the sense of 
geometry (like in Figures ^ and and they most probably should not be interpreted as the excited (i.e., vibrational) 
modes of one of (the segments of) the strings. Such excited modes were studied in the case of a single QCD string in a 
static meson , and in a case of the simply-connected F-shaped string segments in the case of a static baryon [2^ . 
In the sense of such excited mode both energies Eo and Ei correspond to the ground states of the particular string 
configurations visualized schematically in Figures and [21 



III. SOAP-FILM APPROXIMATION 



In the soap-film approximation we evaluate the average of the Wilson loops taking into account only the area 
contributions 

{Q\Wi . . . Wn\0) =^ exp(-a5,), (12) 

i 

where the sum is taken over all topologically distinct and locally minimal surfaces with the area Si which are spanned 
on a given set of oriented contours. The integer pre-factors ki count the number of string configurations which can 
be spanned on the given set of loops. Finally, a is the universal string tension. 

In general, Eq. (|12|) can be justified either in the strong coupling regime (discussed in Ref. for the multi-quark 
interaction) or in the infrared region, where the perimeter-like contributions to the Wilson loops can be neglected. 
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Moreover, the soap-film expansion - applied in the infrared region - coincides with the field-strength correlator 
expansion performed in Ref. [23| to study interactions of the Wilson loops in the Gaussian approximation. 

The essential point is that in Eq. H12|l we neglect the string fluctuations, the contribution of which is small compared 
with the leading area term. Note that in the discussed approximation (|12ll any local string interactions {e.g. via 
exchange of a particle) are neglected. On the other hand, the effect of re-alignment is still encoded in Eq. H12|l. 
therefore this representation ideally suits our needs as a simple theoretical model. 

The expansion of the matrix elements ISJ,® is shown in Table Hll Overall normalization of the matrix elements 
can easily be obtained from the normalization of the basis states. The color factors 1/Nc and in the non-planar 

contributions to the diagonal matrix elements appears in the strong coupling expansion and in the field-strength 
approach [2p|. 

Decomposing the surfaces into the spatial and the time-like faces, 

exp(-a5) = l[ expi-aS,) , (13) 

faces 

and identifying the spatial contributions with the basis states overlap .9, one derives 
where 

fr = exp(— 2(Trt) , (15) 

and 

J. I exp(— 2(7c?t) opposite orientation ^^^^ 

1 exp(— (T(\/3d + r)t) same orientation 

The characteristic equation corresponding to Eq. ||SJ|, 

det(S(t) - AT) = , (17) 

where T is defined in ©, has the following roots 

1 



Ao,l = -(/d + fr± Vifd - fr? + ^g^fdfr) ■ (18) 

Finally, the energy levels are obtained with the help of Eq. (|11|) . 

We are interested in the two regions of the parameter space: (i) the region near the flip-flop point r = d and (ii) 
the "interaction" region corresponding to the two relatively long strings placed near each other, d <C r. 

d — dc At the flip-flop point fr = fd = f ■ Then Ea. (|18|l becomes 

Ao4U=d. =/(l±ff), (19) 
giving the levels splitting (valid for both string orientations) 

{E,-E,)U^,^.^\ln\±Ji^^Ji±^. (20) 
t 1-g t 

Thus the level splitting is just a finite t effect in the case when the basis of the Hamiltonian eigenstates is 
truncated (the true Hamiltonian eigenvalues do not depend on t in the full basis). 

d ^ ?■ At the interaction region the energy splitting can be obtained by an expansion of the square root in Eq. ()18|l in 
powers of fr/ fd, 

Ao = fd{l + fr/fd-g'' + o{fr/fd)} (21a) 

Al = fr{l-9^ - frlfd-g\l-g^)+o{{frlfdf)}. (21b) 
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Then the leading contribution to the energy levels takes the following form (the orientation is explicitly indi- 
cated), 

I j 2(7d opposite , . 

I cr(voa + r) same 

E^U^o = 2ar - = 2ar + ^1±^ , (22b) 

leading just to the "classical" string configurations up to the 1/t corrections. In this case the interaction is 
absent as it was expected. 

An exact (uncxpandcd) expressions for the eigenvalues as functions of the parameters d and r are too bulky we 
skip their functional forms here. However, we present the soap-film results in the graphical form in the next Section 
in order to compare them with the lattice data. 

For the comparison of the lattice data and the soap-film predictions, the latter should be improved to cover the 
small d and r regions in a parameter space, for which the Wilson loops are not saturated by area law. A possible 
way to do this systematically is to make in Eq. H13I) a substitution of single Wilson loop averages for time-like faces 
(which is supposed to be an area-law behavior in the soap-film formalism) by the expression 

expi^aSi) ^ (W,) = exp(-y(r,)i) (23) 

where i runs over all plain time-like faces of size x t, bounded by rectangular Wilson loop and (W) is related to the 
static quark-antiquark potential V{r). In fact we already performed this trick when wc used for the spatial faces the 
overlap g(r, d) instead of the area-law expression which would originate from the soap-film formalism. 
The equations l|15(l . (|16|l now look as: 

fr - exp(-2\/(r)i), (24) 

and 

^ _ I exp{—2V{d)t) opposite orientation ^^g-j 

~ I exp(-(4y(d/V3) + V{r - d/V3))t) same orientation 

All other expressions remain unchanged. In order to perform the improvement we evaluate numerically the static 
potential V{r) and the overlap function g{r, d). Below we use the improved soap-film formulation instead of the naive 
area-law expansion. 



IV. LATTICE DATA 



In order to obtain energy eigenvalues and eigenvectors one needs to know the evolution matrix elements Sab given 
in Eq.QjiO, the overlap of the basis states g(r, d) and the static potential V{r). Since all mentioned observables are 
expressed in terms of the Wilson loops correlators, we perform a lattice computation which is straightforward. The 
details of numerical setup and computation details are given in Appendix 1X1 

The analysis is performed similarly to the analytical case of the soap-film expansion presented in the previous 
Section. The energy levels of oppositely aligned and co-aligned strings are provided as functions of parameters r, d in, 
respectively. Tables IVl and IVII of Appendix IbI 

In Fig. Owe show the energy eigenvalues i?o.i(r, d) as the functions of the string separation d for relatively large 
length of the strings, r = 6a = 0.84 fm. In these figures the energy of the two strings 2V{r) is subtracted to focus 
on the interaction between the strings (this subtraction corresponds to the normalization factor in the denominator 
of Eq. Both figures clearly show the string re- arrangement at some point. We found that the gap between the 

states in the flip-flop points exists in all the studied cases and it is of the same order for the co-aligned and oppositely 
aligned strings. The soap film prediction is shown in Fig. |21by solid lines®. 



Note that the soap-film inputs (the overlap g and the static potential V) were calculated for the on-axis values of parameters, which 
correspond to integer-valued distances in units of the lattice spacing a. In order to achieve a more frequent sampling of the soap-film 
predictions, the values and the statistical errors of g{r,d) and V{r) were interpolated to the non-integer r/a and d/a with the step 
a/10. The ambiguity of such interpolation does not infiuencc the comparison of the soap-film with the data, since the data itself was 
calculated only for on-axis parameters. 
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FIG. 3: The energy eigenvalues Eo,i{r,d) as the functions of the string separation d at the fixed string length r = 0.84 fm. 
The left (right) plot corresponds to the oppositely aligned (co-aligned) strings. The points corresponds to the numerical data 
obtained in simulations while the lines are the soap-film predictions. The statistical errors of both data points are of the order 
of the symbol size. The error in estimation of the soap film prediction is of the same magnitude. Note that in all our numerical 
results the perturbative contributions are suppressed by the HYP procedure 



The interaction energies - shown in Fig. |31 ~ are in a very good agreement with the soap-film prediction discussed 
in previous Section. To compare the soap-film model with the numerical data we plot the difference between the 
predicted and measured energies for the oppositely (likewise) aligned strings of the length r = 0.84 fm in Fig. ^ 
(Fig. |3J). One can clearly see that in the case of the opposite orientations (Fig. the deviation is very small (of 
order of lOMeV) compared to the strings energies (approximately 1 GeV) and is enhanced in the fiip-flop point. The 
deviation is -10 MeV at the flip-fiop point both for the ground state and for most length of the string in the exited 
state. At the lowest distances the deviation is of the order of 50 MeV for the exited state. This might be due to 
overlap of the cores of the strings at these distances. 

Besides the gap in the re-arrangement point, the energy of the ground state of the strings configuration with the 
same orientation deviates from the soap-film prediction at the string separations smaller than dc — 0.4 fm, Fig. |S1 
The exited state has a bit different slope than the soap-film prediction. The deviations are of order of 100 MeV which 
is still at 10% level compared to the energy of a string. The 10% deviations observed in Fig.|31can be explained by 
imperfection of the initial and final states of the string. Indeed, the initial string is made of segments of the Schwinger 
lines, the length of which is of order of the physical string width. 

In Fig. we show the energies Eo,i{r, d) as functions of r for d = dc- One can see that at the flip-fiop point the SF 
model works very nice for all distances r apart from small ones. The discrepancy at small r might be due to effects 
of HYP procedure. Comparison of data obtained with t = 4 and t = 6 Wilson loops show that for large enough r the 
level splitting indeed decreases with increasing t, as predicted by Eg. 1201 From our data we can not conclude whether 
the splitting goes to zero in the limit i — > oo. 

Let us stress that the accuracy of the soap- film expansion is amusingly high taking into account the scales of the axis 
in Figures 21 and [S] and the fact that no fits are made to achieve this nice agreement between the soap-film predictions 
and the numerical results. 

In addition to the energy we study the eigenstates of the Hamiltonian. The energy eigenstates |ft.o,i) are shown as 
the functions of d for fixed r = 0.84 fm in Fig. [7| One can observe the flip-flop behavior in the coefficients hiA and 

hiB, 

\h,) ^ h,A{r,d)\A) + h,B{r,d)\B) , (26) 

which describe the expansion of the ground and excited states in terms of the original states \A) and \B). Note that 
the re-arrangement of the strings is visible simultaneously both for the ground and for the excited states. The absolute 
value of the coefficients may be larger than unity since the states |^) and \B) are not orthogonal to each other. 

The soap-film prediction is also shown in Fig.[7| As in the case of the energy levels, it is almost perfectly consistent 
with the data. 
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FIG. 4: The difference between 5E = E — Eso&p the measured energy levels and the prediction of the soap-film approximation 
for the oppositely aligned strings of the length r = 0.84 fm. The points represent 5E while the band represent the statistical 
errors of soap-film prediction. Note the MeV scale of vertical axis compared to the GeV scale of interaction energy in Fig. |31 
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FIG. 5: The same as in Fig. |l|but for the likewise oriented strings. 
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FIG. 6; Ground and excited energies for oppositely aligned strings at the flip-flop point vs. the distance r between the strings 
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FIG. 7: The same as in Fig. |3 but for the expansion of the energy eigenstates in the basis of the initial states \A} and \B). 
The left (right) plot corresponds to the oppositely (co-aligned) aligned strings while top (bottom) row corresponds to ground 
(excited) states. 



V. INTERPRETATION AND CONCLUSIONS 



We studied the interaction energy of the tetraquark system focusing on the interaction between the string states 
which appear in this quark system. Using a two-state approximation we calculate the eigenf unctions and eigenvalues 
corresponding to the string states of different geometries as the function of the positions of the tetraquark constituents. 
The data is in a good agreement with the prediction of the simple soap-film model which can be considered as an 
analog of the strong coupling expansion or as a zero order (Gaussian) approximation in the field strength correlation 
method 0, . We observe a small deviation from the soap-film prediction at the flip-flop positions of the tetraquark 
constituents where the classical energies of the different string configurations (allowed by the (anti-)quark arrangement) 
coincide. This deviation is very small (of the order of 10 McV for the ground state) compared to the string energies 
(which are typically of the order of 1 GeV). This means that in the realignment point a possible "additional" (i.e., 
is not caused by the flip-flop) interaction between the string segments is very small. Our conclusion is in agreement 
with similar observation made in the SU(2) gauge theory [13j |. 

It is interesting to discuss the relation of our results to the dual superconductor model in QCD which also predicts a 
certain interaction between the QCD strings. The dual superconductor model was invented to describe the confining 
properties of the Yang-Mills vacuum. The main role in the dual superconductor is played by special configurations 
of the gluonic field which are called Abelian monopoles. Such a configuration has a Dirac-like monopole singularity 
in its center after a suitable Abelian gauge is fixed. The importance of the monopoles is stressed by the observation 
of the Abelian dominance implying that the tension of the chromoelectric string is dominated by the Abelian 
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monopole contributions. The model predicts that the QCD strings are interacting via exchange of the dual gauge 
boson and the dual Higgs fields. 

Various lattice simulations indicate that the non-Abelian gauge theories (both SU(2) and SU(3)) correspond to the 
border between type-1 and type-2 dual superconductors. At this border the masses of the both dual fields coincide 
and the parallel strings with the same orientation are not interacting at all while the string with opposite interactions 
are always attracting. This interaction goes beyond the soap-film prediction which works very well as we have just 
noted. Moreover, we see that even at smallest distances the deviation of the ground state energy from the soap-film 
picture is very small, and therefore the traces of the additional inter-string interaction are not found. 

On the other hand it is expected that the masses of both fields are of the order 1 GeV which, in turn, is the 
scale of the lattice spacing used in our model. Since the inverse mass of the dual fields is of the same order as the 
lattice spacing, our data can not either confirm nor reject the validity of the applicability of the dual superconductor 
approach as a possible model for the string-string interactions in QCD. 
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APPENDIX A: NUMERICAL SETUP 

We use a set of quenched SU (3) gauge field configurations generated by MILC from the Gauge Connection (2^ . 
The parameters of these configurations are shown in Table IIIII To enhance statistics the combination of spatial 
APE smearing and time-like HYP smearing was used. The corresponding definitions and parameters are given in 
Appendix lA II for the APE smearing and in Appendix I A 21 for the HYP smearing. 



Action: Tadpole and Symanzik improved 

Coupling: (3 = 8.0 

Size: 20^ x 64 

Number''': full set: 390 

reduced set: 150 
Lattice spacing''': a = 0.14fm 

Smearing parameters'^': asm = 0.55, Nsm = 25 
HYP parameters'*' : ai = 1.8, 02 = 0.6, 0:3 = 0.3 



TABLE III: The simulation parameters. "' The number of uncorrelated configurations. Most data was obtained on the full 
set of configurations, and the data for the co-aligned strings was obtained on a reduced set of configurations; ''' Assuming 
= 440 MeV; '=' See for the definitions; '*' See for the definitions. 

The expectation values for all Wilson loops were measured for the on-axis geometries discussed in the body of the 
text. 



1. APE Smearing 

We use the traditional 3D APE smearing 30] in order to enhance the overlap of the states, created by the gauge 
transporters, with the states of the physical confining strings. It is defined as an iterative fields transformation, one 
iteration has the form: 

<r = ( (1 - c.sn.)U.. + ^ E U.,U^+i.ul,^^ I , (Al) 

"^^^^ ' SU{3) 
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where i,j = 1 . . .3 and (•)s(7(3) rneans a projection onto an 5*17(3) element. Equation IjAip is applied to all spatial 
links and then A^sm iterations are performed. The parameters cksmi are tuned to compromise CPU time with an 
enhancement of the "almost spatial" Wilson loop R x T ~ 10 x 1, see Ref. for details. Used values are: 

a,™ =0.55 A,™ = 25 (A2) 

The overlap of the state, created by the smeared transporter, and the ground state of the confining string ^ is consistent 
with unity {i.e., a complete overlap) within the errors for all values of the string length used. 



2. Hypercubic Blocking 

The hypercubic blocking (HYP) 1^3 was proved to be an efficient tool for the static potential measurements. It is 
defined as a three step transformation of the time-like links: 



U^Up = aiC/.,M+ E U,,nU.+f,,,Ul^i, \ (A3a) 



±p=/^p,i' 



5(7(3) 



C^ifr = «3f/.,.+ E t^i'^;.<Wf^i+L;J ' (A3c) 
V / sum 

where are parameters of the HYP. We do not iterate the procedure in order to preserve locality. To tune the HYP 
cn,2,3 we scan the parameter space for the valley of a maximal enhancement of the "almost time-like" Wilson loop 
T X i? = 10 X 4. The obtained values of the parameters at are: 

= 1.8 a2 = 0.6 = 0.3 (A4) 

Note that these values are quite different from "canonical" ones, obtained by an indirect method, see Ref. 1^3 for 
details. 

The principal effect of the HYP is a suppression of perturbative contributions to the 2- and 4-quark potentials. For 
example, the static potential V{r), calculated after the application of the HYP procedure with parameters chosen 
in this paper, has zero quark self-energy term, a small Coulomb term and a linear term with the unchanged string 
tension. 

The APE and the HYP smearings do not interfere with each other. Therefore the transformations (|A1|) and (|A3|) 
are both applied to the initial links. 



APPENDIX B: 2 AND 4-QUARK POTENTIALS 
1. Remarks on systematical errors 

In our paper we use a rather small time extension t = 4a of the Wilson loops. To check the magnitude of the related 
systematical errors additional measurements were performed for a larger t with a smaller statistics. 

We calculated the 2-quark static potential using the Wilson loops with time extension t' = 5 — 8a and the 4-quark 
potential for the Wilson loops with t' — 6a. 

It occurs that the ^-dependent contribution to V{r) is 0.1% from t-independent contribution, see Table HVl The 
same ratio for the 2-string energy levels is of order of 10^^ (not shown in tables). Thus we estimate a possible 
systematic error coming from finiteness of t to be not higher than 0.1 — 1%. 



Do not confuse this ground state of one string with the ground state of the two strings system, discussed through the paper. 
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Numerical data 



tjr 



10 



aV(:r,t = 4) 

0.057987(9) 0.16662(3) 0.30149(7) 0.4209(1) 0.5331(1) 0.6420(2) 0.7493(3) 0.8560(5) 0.9625(7) 1.069(1) 



0.9969(2) 1.0002(3) 1.0016(3) 1.0018(4) 1.0015(5) 1.0010(7) 1.0004(8) 1.000(1) 1.000(1) 1.000(2) 

0.9947(2) 1.0002(3) 1.0026(3) 1.0031(5) 1.0026(6) 1.0017(8) 1.001(1) 0.999(1) 0.999(2) 1.000(3) 

0.9932(2) 1.0003(3) 1.0034(4) 1.0039(5) 1.0032(7) 1.002(1) 1.001(1) 1.001(2) 1.003(4) 1.009(8) 

0.9920(2) 1.0003(3) 1.0038(4) 1.0042(6) 1.0034(9) 1.002(1) 1.000(2) 1.000(4) 1.001(8) 1.00(1) 



TABLE IV: Static quark-antiquark potential Vij) (in lattice units) and its t-dependence 



d/r 



10 



2 
3 
4 
5 
6 
7 
8 
9 
10 



0.40692(3) 
0.46497(4) 
0.46887(4) 
0.46951(2) 
0.46964(3) 
0.46967(3) 
0.46968(9) 
0.46968(3) 
0.46969(2) 



0.46497(4) 
0.78368(8) 
0.84370(8) 
0.84875(5) 
0.84960(6) 
0.84979(6) 
0.84984(6) 
0.84984(6) 
0.84990(4) 



0.46887(3) 

0.84369(7) 

1.1268(1) 

1.18117(9) 

1.1854(1) 

1.1862(1) 

1.1864(1) 

1.1864(1) 

1.18658(6) 



Ground 

0.46951(2) 

0.84875(5) 

1.18117(9) 

1.4540(1) 

1.4992(2) 

1.5021(2) 

1.5026(2) 

1.5027(2) 

1.5029(1) 



state Eii(r, 

0.46963(3) 

0.84960(6) 

1.1855(1) 

1.4992(2) 

1.7743(3) 

1.8081(3) 

1.8096(3) 

1.8097(3) 

1.8100(2) 



d), GeV 

0.46967(3) 

0.84979(6) 

1.1862(1) 

1.5021(2) 

1.8082(3) 

2.0910(6) 

2.1122(5) 

2.1120(5) 

2.1129(4) 



0.46968(3) 

0.84984(6) 

1.1864(1) 

1.5027(2) 

1.8096(3) 

2.1122(5) 

2.403(1) 

2.411(1) 

2.4142(8) 



0.46968(3) 

0.84984(6) 

1.1864(1) 

1.5027(2) 

1.8097(3) 

2.1120(5) 

2.412(1) 

2.706(3) 

2.715(2) 



0.46969(2) 

0.84990(4) 

1.18657(6) 

1.5029(1) 

1.8100(2) 

2.1129(4) 

2.4142(8) 

2.715(2) 

3.003(3) 



2 
3 
4 
5 
6 
7 
8 
9 
10 



0.58197(6) 

0.9068(1) 

1.2421(3) 

1.5558(5) 

1.851(2) 

2.124(3) 

2.366(6) 

2.56(1) 

2.71(1) 



0.9069(1) 

0.9561(1) 

1.2273(2) 

1.5330(3) 

1.829(1) 

2.111(2) 

2.375(5) 

2.613(9) 

2.82(1) 



1.2424(3) 

1.2274(2) 

1.2711(2) 

1.5264(2) 

1.8218(7) 

2.113(1) 

2.397(3) 

2.671(7) 

2.929(9) 



Exited 
1.5563(5) 

1.5332(3) 

1.5264(2) 

1.5638(2) 

1.8207(4) 

2.1153(9) 

2.408(2) 

2.699(5) 

2.982(7) 



state E\{r 
1.852(3) 
1.829(1) 
1.8216(7) 
1.8205(4) 
1.8500(4) 
2.1168(7) 
2.414(2) 
2.709(4) 
3.005(6) 



d), GeV 
2.127(3) 
2.111(2) 
2.113(1) 
2.1147(9) 
2.1165(7) 
2.1367(6) 
2.416(1) 
2.715(3) 
3.017(5) 



2.370(6) 

2.376(4) 
2.397(3) 
2.408(2) 
2.413(2) 
2.416(1) 
2.424(1) 
2.712(3) 
3.014(5) 



2.57(1) 

2.616(9) 
2.671(7) 
2.699(5) 
2.710(4) 
2.715(3) 
2.713(3) 
2.715(3) 
3.009(5) 



2.72(1) 

2.83(1) 

2.932(9) 

2.982(7) 

3.004(6) 

3.017(5) 

3.014(5) 

3.010(5) 

3.030(4) 



TABLE V: Energy levels of oppositely aligned strings as functions of geometry r, d (in lattice units, a = 0.14 ftn). 
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d/r 


4 


6 


8 


10 


4 


6 


8 


10 




Ground state Eo{r,d), GeV 


Exited state Ei{r,d), GeV 


2 


0.8829(3) 


1.1943(4) 


1.4950(7) 


1.795(1) 


1.347(1) 


1.983(4) 


2.58(1) 


3.12(5) 


3 


1.0909(3) 


1.4450(5) 


1.7509(7) 


2.052(1) 


1.3689(7) 


1.926(2) 


2.495(6) 


3.03(2) 


4 


1.1683(3) 


1.6571(5) 


1.9856(9) 


2.291(2) 


1.5284(8) 


1.923(1) 


2.470(3) 


3.03(1) 


5 


1.1828(2) 


1.7796(5) 


2.226(1) 


2.549(2) 


1.764(1) 


2.038(1) 


2.469(2) 


3.028(9) 


6 


1.1855(2) 


1.8036(4) 


2.382(2) 


2.800(5) 


2.002(2) 


2.251(2) 


2.554(3) 


3.052(9) 


7 


1.1862(1) 


1.8085(4) 


2.411(2) 


2.995(7) 


2.234(3) 


2.500(3) 


2.789(4) 


3.14(1) 


8 


1.1864(1) 


1.8097(4) 


2.417(1) 


3.035(8) 


2.445(4) 


2.746(5) 


3.047(9) 


3.36(2) 


9 


1.1864(1) 


1.8100(3) 


2.421(1) 


3.042(8) 


2.641(7) 


3.00(1) 


3.32(2) 


3.60(4) 


10 


1.18655(8) 1.8097(3) 


2.416(1) 


3.033(7) 


2.814(8) 


3.24(2) 


3.57(3) 


3.82(7) 



TABLE VI; The same as in Table IVI but for the co-aligned strings. 
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